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Abstract 

Primordial non-Gaussianity has emerged as one of the most promising probes of the inflationary 
epoch. While the cosmic microwave background and large-scale halo bias currently provide the 
most stringent constraints on the non-Gaussian parameter Jnl, the abundance of dark matter 
halos is a complementary probe which may allow tests of Gaussianity which are independent of 
the precise form of non-Gaussian initial conditions. We study the halo mass function in A^-body 
simulations with a range of non-Gaussian initial conditions. In addition to the usual fjsiL model, we 
consider g]\fL^^-type non-Gaussianity and models where the 4-point amplitude t^l is an independent 
parameter. We introduce a new analytic form for the halo mass function in the presence of primordial 
non-Gaussianity, the "log-Edgeworth" mass function, and find good agreement with the A^-body 
simulations. The log-Edgeworth mass function introduces no free parameters and can be constructed 
from first principles for any model of primordial non-Gaussianity. 



1 Introduction 



One of the most exciting prospects in observational cosmology is the opportunity to constrain the 
physics of inflation, thereby probing energy scales which are far beyond the reach of accelerator 
experiments [1, 2, 3, 4, 5, 6, 7]. Current cosmic microwave background (CMB) data provides strong 
evidence for a spatially flat universe with small density perturbations drawn from a nearly scale- 
invariant power spectrum, in accord with inflationary predictions [8]. Nevertheless, distinguishing 
between microphysical models on the basis of the scalar power spectrum alone remains a challenge. 
For single-field, slow-roll inflation, higher-order non-Gaussian statistics of the curvature perturbation 
are unobservably small [9, 10, 11]. However, there are broad classes of inflationary models - those 
that violate slow-roll, have multiple fields, or modified kinetic terms for example - that can generate 
observable levels of non-Gaussianity (see for example, [12, 13, 14, 15] and references therein). A 
detection of non-Gaussianity would therefore rule out single-field, slow-roll inflation and could also 
be a powerful discriminator between these alternative scenarios. 

At present the tightest constraints on non-Gaussianity in $(x), the primordial curvature per- 
turbation^, come from constraining the amplitude of several higher-point "shapes" inspired by dif- 
ferent inflationary scenarios [16, 17, 18, 19, 20, 21, 22]. For instance, in the so-called local model 
[23, 24, 25, 26, 27], the initial curvature is a non-Gaussian held deflned through 

$(x) = $g(x) + Inl i^ci^f - i'^D) + 9NL i'^ci^f - 3(<1>|)<I>g(x)) + . . . (1) 

where $g is a Gaussian field and fNL,gNL are free parameters. The WMAP constraints on these 
parameters are -10 < /nl < 74 [8] and -7.4 x 10^ < ^atl < 8.2 x 10^ [28] or -12.34 x 10^ < qnl < 
15.58 X 10^ [29] at 95% confidence. This model also generates a scale-dependent signature in the 
bias of dark matter halos [30, 31, 32, 33] that allows for competitive constraints from low-redshift 
data: -29 < Jnl < 70 [34], and -3.5 x 10^ < qnl < 8.2 x 10^ [35] at 95% CL. The Planck CMB 
satellite is expected to achieve 1-a errors that are smaller by a factor of 3-5 [25, 36]. 

Signatures of primordial non-Gaussianity can also appear the abundance of dark matter halos 
[37, 38, 39, 40, 41, 42, 43, 44]. For instance, positive (negative) skewness in the density field will 
tend to increase (decrease) the number of very high mass halos hosting galaxy clusters. The number 
density of halos as inferred from the cluster mass function has been shown to be a probe of primordial 
non-Gaussianity which is complementary to the CMB [45, 46, 47, 48, 49, 50, 51] (see also e.g. 
[52, 53, 54, 55] for investigations of weak lensing as a probe of primordial non-Gaussianity). The 
mass function is sensitive to cumulants beyond the 3-point function but relatively insensitive to the 
precise shape of the A^-point functions. Therefore, the mass function can constrain non-Gaussianity 
without prior knowledge of template shapes, but is less powerful for discriminating between forms 
of non-Gaussian initial conditions. At present evidence of any primordial non-Gaussianity would be 
extraordinary, so it is useful to obtain observational constraints from a variety of methods. Recently 
there have been hints of an overabundance of high-z massive clusters [56, 57, 58, 59] which can be 

^In Eq. (1) and throughout the rest of this paper, we have defined a curvature $ = |(^, wliere is the primordial 
curvature fluctuation conserved on super-horizon scales. This notation is conventional in studies of primordial non- 
Gaussianity, although possibly confusing since the Bardeen curvature $h is equal to |^ at early times when Eq. (1) 
is applied but |(^ in the matter dominated era, long after the non-Gaussianity in Eq. (1) is imprinted. 
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interpreted as evidence for primordial non-Gaussianity (but note that [60, 61, 62] find consistency 
with a Gaussian mass function). This picture will undoubtedly sharpen in the near future with 
improved mass function constraints from experiments such as SDSS [63], Chandra [64], ACT [65], 
SPT [66], and Planck [67]. 

In this paper we study the halo mass function in A^-body simulations with non-Gaussian initial 
conditions with two forms of non-negligible trispectra: (i) initial conditions where the connected 
4-point function has the same "shape" as that from the f'^^ terms in Eq. (1) but boosted in am- 
plitude relative to the 3-point function (the ^^tnl" model ) and (ii) initial conditions with the gpfL 
contribution from Eq. (1) much larger than that from /nl- This extends current high resolution 
studies of the halo mass function with pure /AfL-type initial conditions [30, 68, 69, 70, 71]. See also 
[35] for studies of the mass function from A-body simulations with a cubic term g^vL in the initial 
conditions, and [72] for A-body simulations with initial conditions with more general primordial 
bispectra. 

We propose an analytic halo mass function, the "log-Edgeworth" mass function, which accurately 
describes our simulations for a wide range of /jvl, gNL, and tnl- The log-Edgeworth mass function 
is conceptually similar to the second-order Edgeworth mass function from [46] , but is a better fit to 
A-body simulations in cases where the two disagree, in particular for the high-mass limit, where the 
Edgeworth mass function breaks down. 

Throughout this paper we use the WMAP5+BA0+SN fiducial cosmology [73]: baryon density 
flbh'^ = 0.0226, cold dark matter (CDM) density ^ch^ = 0.114, Hubble parameter h = 0.70, spec- 
tral index rig = 0.961, optical depth r = 0.080, and power-law initial curvature power spectrum 
PP^(fc)/27r2 = A2(A:/ikpiv)"^"^ where = 2.42 x 10"^ and fcpiv = 0.002 Mpc"!. 

In §2 we introduce the generalized local non-Gaussian initial conditions considered in this paper. 
In §3 we discuss prescriptions for analytic mass functions that describe the effects of non-Gaussianity 
through the cumulants of the density field smoothed on scale M, and also present fitting formulae for 
the smoothed skewness and kurtosis. The A-body simulations and a comparison with analytic mass 
functions is presented in §4. Concluding remarks are given in §5. Appendix A contains a discussion 
of calculations of the smoothed skewness, kurtosis and the 0{ffj^) correction to the variance. 

2 Non-Gaussian Initial Conditions 

The simplest model of primordial non-Gaussianity is the "local" type, in which the initial curvature 
$ is given by 

<I>(x) = ci>G(x) + /;VL(<&G(x)2-(<l>|)) (2) 
and the 3-point and connected 4-point functions of the initial curvature are given by 

(«>(ki)$(k2)$(k3)) = fNLP^{ki)P^{k2){27rf6D (^k,) + (5 perm.) + 0(7^^) (3) 

(«>(ki)$(k2)«'(k3)$(k4))e = 2f^Mki)P'S>{k2)P^{\ki+k3\){2TTfSD (j^ki) + (23 perm.)+0(/^i) 

(4) 

where do is the Dirac delta function and P$ is the power spectrum of ^g- This type of non- 
Gaussianity is generated, e.g. in the curvaton model, in which there is a second light field present 
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during inflation (the curvaton) that decays after inflation has ended and generates the primordial 
curvature perturbation [74, 75, 76, 77, 78]. We assume throughout that the curvaton decays before 
dark matter freezout, so that no dark matter isocurvature mode is generated. 

In this section, we will describe two sets of non-Gaussian initial conditions which generalize 
Eq. (2) and give rise to a trispectrum in the density field that is large enough to affect the halo mass 
function. 

2.1 Equal Power from the Curvaton and Inflaton: An Example of t^l [i/nlY 

Most studies of the curvaton model have focused on the case where the curvaton completely dom- 
inates the primordial curvature perturbation. An alternative set-up, considered recently by [79, 
80, 81, 82], supposes that the curvaton and inflaton contribute equally to the primordial curvature 
perturbation. If the inflaton contribution is assumed to be Gaussian and the curvaton contribution 
is non-Gaussian, the statistics of the initial conditions are changed in important ways. In particular 
the higher A^-point functions are boosted in amplitude relative to the bispectrum, and the boost 
depends on the ratio of inflaton-to-curvaton contributions to the curvature. 
More precisely, the primordial curvature in this model is given by 

«>(x) = ^>i(x) + $,(x) (5) 

where and ^>c denote inflaton and curvaton contributions. We assume that and <I>c are uncor- 
related fields with proportional power spectra, i.e. P<^i{k) = -j^r^P^ik) and P^^{k) = j^^P^{k), 
where ^ is a free parameter which represents the ratio of inflaton to curvaton contributions. 
We take to be Gaussian and <I>c to be a field with local non-Gaussianity: 

<I>,(x) = c^,,g(x) + /jvL ($?,g(x) - (^c' g(x))) (6) 

where $c,G is a Gaussian field. 

In this model the power spectrum, bispectrum, and trispectrum of the initial curvature are given 

by: 

(cl.(ki)a>(k2)) = P$(A;i)(27r)35z5(ki+k2) + 0(/^i) (7) 
($(ki)$(k2)$(k3)) = fNL{27Tf6D {Y.^^} ^* (^1 )^<i> (^^2) + (5 perm.) + 0(/^i) (8) 

($(ki)a>(k2)1>(k3)$(k4))e = ^^2(2^)3^^ (J2 k^) P^ (ki) P^ {k2) P,S> + kgl) 

+ (23perm.) + 0(/^^) (9) 



where we have defined 



TNL '- 
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TNL = {^] (1 + a. (10) 



The bispectrum and trispectrum in this model have the same shapes as in the curvaton-dominated 
model considered previously (Eq. (2)), but the coefficients Jnl, tnl are independent parameters. 
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The curvaton-dominated model corresponds to the special case t^l = (f/Afi)^ and the current 
bounds are —6000 < tjs^l < 33,000 at 95% confidence [28]. The factor (6/5) here is conventional 
and has been introduced for consistency with the literature. Our perspective is that primordial 
non-Gaussianity is most conveniently parameterized by the coefficients of the primordial A^-point 
functions, and the specific two-field model in this section is just a mechanism for generating local 
non-Gaussianity with prescribed /nl and tnl- For this reason, we will use /nl and tnl as the 
basic parameters of the model, and treat fi\fL and as derived parameters. 

2.2 Local Initial Conditions with Kurtosis but No Skewness: gj^i 

Another variation on initial conditions with local non-Gaussianity is to consider a case where the 
quadratic term in Eq. (1) vanishes but the cubic term is included [26]. Initial conditions of this 
form can be generated in a curvaton model where the potential for the curvaton has terms that 
are not quadratic and cancellations from these terms set /atl ~ while generating a large g^L 
[83, 77, 27, 84] 

<I>(x) = «I>g(x) + g^L ($|(x) - 3(cD|)c^g(x)) . (11) 

We have found it convenient to include the — 3(<l?g)<I>G(x) term in the definition of the g^L model 
so that the power spectrum of $ will be unchanged to first order in gNL-'^ Note that the expectation 
value {^q) is infrared divergent, but converges in a finite volume. We compute it as a discrete 
sum over Fourier modes k 7^ in our simulation, for consistency with the way the Gaussian initial 
conditions are generated. 

With initial conditions given by (11), the power spectrum and connected trispectrum are given 

by: 

(a>(ki)c|>(k2)) = {2Trf5D{]^i + ]^2)P^{ki) + 0{g%L) (12) 

($(ki)$(k2)$(k3)$(k4))conn. = gNL{27TfSD ( J] k,) P$ (fci )P$ (^2)^-!. (A^s) 

+ (23perm.) + 0(g2^i) (13) 

and all odd A^-point functions vanish. 

3 Halo Mass Functions - Theory 

There are a number of non-Gaussian mass functions in the literature. Most of them are derived from 
some expansion of a non-Gaussian probability distribution function (PDF) for the mass density field 
and either follow the Press-Schechter ansatz [85, 39, 86, 44, 46, 87], or more formally use excursion 
set theory [88, 89, 90, 91]. These mass functions are typically specified by the PDF or iV-point 
correlation functions of the linear density field (or equivalently the initial curvature). On the other 

■^An alternate convention (e.g. [35]) omits this term from the definition; in this case the "bare" power spectrum 
amphtude used to set up the Gaussian field $g will differ from the observed value of that would be inferred 
by measuring observable power spectra. The two definitions will be equivalent in any analysis which marginalizes A^, 
but the first requires less bookkeeping to keep track of the difference between bare and observable power spectrum 
amplitudes. 
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hand, Pillepich, Porciani and Hahn [70] present a fitting formula for the mass function derived from 
non-Gaussian simulations which (in addition to the usual cosmological parameters) depends on /nl, 
assuming qnl = ^ and t^l = (f /a^l)^- These mass functions are all in relatively good agreement 
with each other and with simulations with fNL-type (i.e. Eq. (2)) non-Gaussian initial conditions 
[70, 71]. In this paper, our emphasis will be on analytic mass functions that depend only on the 
cumulants of the variable Sm, the linearly-evolved density fluctuation smoothed on mass scale M 
[46]. However, see [92, 91, 90] for extensions based on excursion set theory which include additional 
terms from so-called unequal time correlators. 

In this section we describe a general formalism for deriving non-Gaussian mass functions, and 
apply it to the case of local non-Gaussianity parameterized by /nl, dNL, and tnl- Our approach is 
conceptually similar to the Edgeworth approach from [46] but differs in some details which we now 
explain. We'll follow the Press-Schechter model [85] which states that the fraction F{M) of volume 
collapsed to objects of mass > M is equal to the probability for 6m to exceed the collapse threshold 
Sc ~ 1.42.'^ In the presence of primordial non-Gaussianity, the 1-point PDF p{Sm) is perturbed, and 
this leads to a change in the mass function which is computable in the Press-Schechter model. 

The Edgeworth expansion is a representation of a general PDF as a power series in the higher 
cumulants of the distribution. Since higher cumulants parameterize deviations from Gaussianity, the 
Edgeworth expansion is most useful in the regime of weak non-Gaussianity, where the series converges 
rapidly. In [46], the Edgeworth expansion for p{6m) was truncated to obtain an estimate of the /nl 
dependence of p{Sm), from which the /nl dependence of the mass function can be calculated. Here, 
we will find it convenient to truncate the series expansion for the quantity ln(F(M)) (rather than the 
quantity p{5m))\ we call the non-Gaussian mass function obtained in this way the "log-Edgeworth" 
mass function. We will find that the log-Edgeworth mass function is a better fit to simulations 
than the Edgeworth expansion for parameter values where the two disagree, in particular for the 
high-mass limit, where the Edgeworth expansion breaks down. 

Another more minor detail is that we keep a few more terms in the series expansion than are 
usually quoted from [46]. We do this for the following reason: the simulations in §4 will show that 
the mass function is TATL-dependent, in the case where tnl is varied at fixed Jnl (and with qnl = 0). 
Therefore, to allow tml dependence, we will keep the term in the Edgeworth series which is first 
order in ttvl. Since ttvl is the same order as /^j^, we will also keep terms of order in order to 
truncate the series consistently. Similarly, we keep terms of first order in g^L but not second order, 
since qnl has the same order as tnl when it appears in A^-point correlation functions. 



3.1 Cumulants 

The linear density field smoothed on scale M is given by 

hi{z) = I ^WM(.k)a{k,zMk) (14) 

''Throughout this paper we will include a correction to the spherical collapse threshold Sc — 1.686, 5c — > Sc^ with 
q = I/V2 which has been shown to give better agreements with simulations [70, 71]. 
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where 

3sin(fci?(M)) 3cos(fcfi(M)) 

^^'^(^^ = (fci?(M))3 (fci;(M))^ 

is the Fourier transform of a tophat window whose comoving radius R{M) = {3M/47rpm)^/'^ encloses 
mass M, and 

, 2k'^T(k)D(z) , , 

is defined so that (5(k, z) = a(A;,2;)$(k) in linear theory. Here, T{k) denotes the matter transfer 
function, and D{z) is the linear growth function normalized to D{z) = 1/(1 + z) at high z. 
The reduced cumulants are defined as 



i^n{M) = ' for > 3 (17) 



where (J^(M) = Note that Kj\f{M) is independent of z as implied by the notation. 

For the non-Gaussian mass function in the next subsection, we will need to know the cumulants 
At3(M) and K4^{M) to leading order in the non-Gaussianity parameters fwL, dNL and tnl- These 
cumulants can be calculated either by Monte Carlo, or analytically by integrating the A'"-point 
correlation function (Eqs. (8), (9), (13)) over the wavenumbers kj with appropriate weighting. The 
reduced cumulants (with the (t(M)^ denominator) are slowly varying functions of M and relatively 
insensitive to the assumed cosmological parameters but slightly cumbersome to compute. Details 
of the calculation are given in Appendix A. Here, we simply quote the results: for k^, we find the 
following to be a good fitting function 

«3(M) « fNL (6.6 X lO-'^) (l - 0.016 In f rj^) ) (18) 



(see also [87]). The cumulant K4 is more subtle: it contains terms proportional to ^tvl and ttvl, and 
the C(t7vl) term formally diverges as the volume of the simulation box is taken to infinity. This 
divergence is well-known and we discuss it in detail in Appendix A. For now, we just quote a fitting 
function for K4 with explicit dependence on the box size L:^ 

At4(M) « (1.6 X 10-^) (^1-0.021 In (^^^3^^^ 



(19) 



where Lq = 1600 h ^ Mpc and A|, = is equal to (8.72 x 10 for the fiducial cosmology from 
§1- 

■'in Eqs. (19) and (21), the factor A| ln(L/Lo) is actually an approximation to the exact infrared divergent behavior 
which assumes (1 — ns)\n{L/Lo) ^1. If this condition is not satisifed then a more accurate approximation can be 
obtained by making the replacement 



A|ln 



fe=4.67/L 
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We will also need the leading non-Gaussian contribution to the variance cr(M)^, which has the 
same order^ (i.e. quadratic in f^L in the case tnl = (S/S/atl)^) as the second term in K4. We write 
the variance as a sum of a Gaussian term (7^(M) and a non-Gaussian term (7^(M)K2(Af ), where K2 
represents the fractional correction due to primordial non-Gaussianity 

C7^(M) = a^(M)(l + K2(M)) (20) 

and find the following fitting function for K2{M) 



K2{M) ~ 



(4.0 X 10-8) ^1 - 0.021 In (^-^^) ) + 4A| In ^ ^ 



(6/5)4/2,^ ^ V \h-^Mr^)) * V^o 



(21) 



where Lq = 1600 h""^ Mpc. Note that K2(M) is also infrared divergent. 

The correct choice of L in the infrared divergent parts of the above cumulants depends on the 
context in which a non-Gaussian mass function is desired. When comparing to A^-body results in §4, 
we will choose L to be the side length of the (periodic) simulation volume. When comparing to the 
observed mass function from real data, a reasonable choice would be L ~ 2R where R ~ 14000 Mpc 
is the comoving causal horizon. The need to choose a value of L may seem unfortunate, but it is a 
generic feature of models with local non-Gaussianity that has nothing to do with mass functions per 
se: such models only make sense when regulated in a finite volume, and physical quantities (such as 
the non-Gaussian contribution to the power spectrum) depend weakly on this choice. 

In Fig. 1, we plot the reduced cumulants K2, K3, obtained by Monte Carlo, with the fitting 
functions shown for comparison. 

3.2 Mass Function Derivation 

In the Press-Schechter model, the fraction of volume F{M) collapsed to halos of mass M is given by 

F{M) = / dup{u,M) (22) 

Juc.{M) 

where p(iy,M) is the 1-point PDF of the variable u = 6M/cr{M), and Uc{M) = 6c/(y{M) with 
5c = 1.42. The mass function n(M) is then given by 

n(M) = -2^F'{M) (23) 

where primes denote derivatives with respect to M. 

In a Gaussian cosmology, we havep(i/,M) = (27r)-i/2 exp(-zyV2)- In a non-Gaussian cosmology, 
we can expand the PDF as a series (the Edgeworth expansion): 

p{^, M) = + M) +p2{u,M) + ...) (24) 



^Strictly speaking, the amplitude of K2 need not be proportional to t^l/ /nl but may be a free parameter. The 
form of K2 given in Eq. (21) assumes initial conditions with /ml and tnl which are implemented using the two-field 
model from 52.1. 
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1012 IQ13 10'* lO's lO's 
M (M^h) 

Figure 1: The reduced cumulants of the density field smoothed on scale M. The points are the 
cumulants measured from Monte-Carlo realizations of the initial conditions, the lines are the fits 
from Eqs. (18)-(21). Note that the smoothed kurtosis has roughly the same M dependence for 
the qnl (dashed, red line) and tnl (solid, blue line) models despite the different "shapes" of their 
trispectra. For K2 and k^^tnl have assume L = 1600/i~^Mpc. 



where 

Pi{u,M) = ]:Ks{M)H3{u) (25) 
o 

P2{i^,M) = ]^K2{M)H2{v) + ^^K^{M)H^{u) + ^K^{MfH^{u) (26) 

and Hn{i') are the Hermite polynomials defined by Hn{u) = (— l)"e'^^/^^|^e~^^/^. Note that pi 
represents contributions which are first order in fi\[L, and p2 represents contributions which are either 
second order in fiy^ or first order in gjy^, tnl- It is worth noting that the explicit mass dependence 
of pi, P2, arising from the mass dependence of the k(M)'s (as well as the box dependence in K2{M, L), 
At4(M, L)) breaks the universality of the mass function. However as can be seen in Eqs. (19), (21) and 
Fig. 1, the box dependent contributions are small and the mass dependence is slight, so universality 
is only weakly broken. 

Plugging into Eqs. (22), (23), we can obtain series expansions for F{M) or n{M). The series 
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F{M) = Fo(M) + Fi(M) + F2{M) is given by 

V2 J 



Fo{M) 
Fi{M) 
F2{M) 



— erfc 
2 



1 



(27r)i/2 
1 

(27r)i/2' 



6 

K2{M) 



H2{uc{M)) 
Hi{v,{M)) + 



K4(M) 



K3(M)^ 



24 ^ " 72 

It will also be convenient to have expressions for the derivatives with respect to M 



F[{M) 
Fi{M) 



(27r)V2 



-iyc{Mf/2 



F^{M) 

4(M) 
21/^ 



6 &v'. 



24 



Hi{uc{M)) 



24i/' 



72 

K3(M)4(M) 



^5(z^c(M)) 



(27) 
(28) 

(29) 

(30) 
(31) 

(32) 



The "Edgeworth" mass function from [46] is defined by truncating the series for F{M) (or for n(Af )). 
At second order in /nl (first order in ^tvl, tj^l) the resulting mass function is 



riNG 



no 



Edgeworth 



1 + 



F[{M) F^{M) \ 
F^{M)^ Fl,{M)) ■ 



(33) 



This mass function has been shown to be in good agreement with simulations for t^l = i^fNi)'^ 
[70, 71]. However, it has a couple of slightly annoying properties that arise at large [/atlI; if truncated 
at first order in fj^L it is unbounded from below as /nl approaches large negative values, if truncated 
at second order the mass function is non-monatonic at the high-mass end for negative Jnl-^ 

Inspired by these issues, we propose a slightly different mass function. If we truncate the series 
for ln(F(M)) rather than F{M): 

(34) 



ln(F(M)).lnFo(M) + ^ + ^-i 



Fi{M) \ 
Fo{M)) 



then we obtain the mass function 
nNG 



riG 



log-Edgcworth 



exp 



X 



Fi{M) F2{M) 1 
Fo(M) + Fo(M) " 2 



Fi{M) V 
Fo{M)) 



(35) 



Fi{M) + F^{M) _ Fi{M)F[{M) _ F^jM) + F2{M) Fi{Mf 

^ 1-!/ / Ti ,r\ T-1 / -71 r\ T-il / n r\ i-t / -n ,r\ 



Fo{M)F^{M) 



Fo{M) 



^We have experimented with using a mass-dependent barrier Uc{M) in Eqs. (27)-(29), chosen so that the Gaussian 
mass function derived in this way agrees with the Sheth-Tormen or Warren mass functions [93, 94] but found it did 
not improve the Edgeworth mass function. 
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We will refer to this non-Gaussian mass function as the "log-Edgeworth" mass function. 

The Edgeworth (33) and log-Edgeworth (35) mass functions agree to second order in /at^ (or 
first order in qnl, tnl) but differ in the way the low-order derivatives are extrapolated to finite 
values oi /nl, Qnl and tnl- As we will see in the next section, the log-Edgeworth mass function is 
a better fit to A^-body simulations in the regime where the two disagree, and also provides sensible 
asymptotic behavior in the limit of large halo mass. In principle, if we kept all terms in each series 
expansion the Edgeworth and log-Edgeworth mass functions would agree. It is interesting to note, 
that at lowest order \n f^L-, dNL-, the ratio nNc/f^G as predicted by the log-Edgeworth mass function 
reduces to that from the Edgeworth mass function in the limit Uc{M) << 1 and the "MVJ" [44] 
mass function for fc(M) >> 1.^ 

Note that we have written both the Edgeworth and log-Edgeworth mass functions as expres- 
sions for riMG/f^Gi to separate the issue of analytically describing the Gaussian mass function no 
(see e.g. [95, 94, 96]) from the issue of describing the fractional correction due to primordial non- 
Gaussianity. 

In this section, we have presented the log-Edgeworth mass function in maximum generality in 
order to provide a framework which can be adapted to models of non-Gaussianity not considered 
in this paper. If all that is desired is to compute the mass function for specified values of Jnl-, 
qnl and ttvl, this can be done straightforwardly by using the fitting functions in Eqs. (18)-(21) to 
compute Kj(M), then using Eqs. (27)-(32) to compute Fi{M) and Fl{M), and finally using Eq. (35) 
to compute nNc/nc- 

4 Halo Mass Function from A^-body Simulations 

To study the mass function, we performed collisionless A^-body simulations using the GADGET-2 
TreePM code [97]. Simulations were done using periodic box size Rhox = 1600 h^^ Mpc, particle 

1 /3 

count Np = 1024'^, and force softening length Rg = 0.05(i?box/Ap ). With these parameters and 
the fiducial cosmology from §1, the particle mass is nip = 2.92 x 10^^ Mq. 

We make non-Gaussian simulations of the initial curvature ^ in two cases: either (1) taking 
9nl = and nonzero /nl, tnl, or (2) taking /nl = tnl = and nonzero qnl- In the first case, we 
use the two-field model from §2.1. Starting from input parameters fiy^ and tnl, we obtain f^^ and 
^ using Eq. (10), then simulate Gaussian fields $j and ^c,G with power spectra ^^^^ and j^^^P^ 
respectively. The total initial curvature is then given by <I> = <I>j + <I>c, where <I>c is given by Eq. (6). 
In the second case (/atl = ttvl = with qnl / 0), we simulate a non-Gaussian $ using Eq. (11). 
We do not generate initial conditions with nonzero values for all three parameters /nl, 9nl, tnl, 
although it would be easy to extend the two-field model to apply in this generality. 

Given a realization of the initial curvature initial particle positions and velocities are generated 
as follows. First, we apply the transfer function T{k), computed using CAMB [98], to obtain the 
Newtonian potential at the initial redshift Zjni = 100 of the simulations. Then we obtain initial 
particle positions using the Zeldovich approximation [99]. (At Zini = 100, transient effects due to 
use of this approximation should be negligible [100].) 

^We are grateful to Vincent Desjacques for pointing this out. 
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Figure 2: Comparison of the Edgeworth (Eq. (33)) and log-Edgeworth (Eq. (35)) mass functions for 
non-Gaussian initial conditions with nonzero Jnl and tnl- For tnl = {^Inl)"^ (i-e. perturbations 
generated entirely by the curvaton) they both provide reasonably good fits. For tnl = '^{^Inl)'^ 
(i.e. equal power from the curvaton and inflaton) the log-Edgeworth mass function is in better 
agreement. 
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Figure 3: Comparison of mass fmictions with tnl = {^fNi)"^ and tnl = 2(|/ArL)^ for: {a)fNL = 
+250 (b)/ArL = +500 {c)fNL = -250 and (d) Jml = -500. This plot is intended to isolate the effect 
of the primordial trispectrum on the halo mass function, since the two models being compared have 
equal bispectra but different trispectra. The curves are the "Edgeworth" mass function in Eq. (33) 
and the "log-Edgeworth" mass function presented in Eq. (35). 
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Figure 4: Comparison of the Edgeworth (Eq. (33)) and log-Edgeworth (Eq. (35)) mass functions for 
initial conditions with a qnl^^ term (Eq. (11)) (here /nl = tnl = 0). 



After running the A^-body simulation, we group particles into halos using an MPI parallelized 

— 1/3 

implementation of the friends-of-friends algorithm [101] with link length Lfof = 0.2i?boxAp . For 
a halo containing NpoF particles, we assign a halo position given by the mean of the individual 
particle positions, and a halo mass given by 

rrih = nup (iVpoF - A^fof) > (36) 

the second term is recommended in [94] to minimize particle resolution artifacts in the mass function.® 
In Fig. 2 we study the mass function n[M) for a few choices of {/nLiTnl). It is seen that the 
mass function mainly depends on /tvl, but some TATL-dependence can also be seen, particularly at 
high redshift. Recall from Eq. (8) that the 3-point function is independent of ttvl, so the ttvl- 
dependence at high redshift indicates that the skewness of the density field (in addition to the linear 
power spectrum) is insufficient to describe the non-Gaussian mass function. However, the four-point 
function (Eq. (9)) and the 0{tnl) correction to the power spectrum (Eq. (21)) do depend on tnl, 
so our simulations indicate that these higher cumulants are relevant for determining the abundance 
of high-mass halos. 

The Edgeworth and log-Edgeworth mass functions given in Eq. (33) and Eq. (35) include the 
first TATi-dependent terms in the series expansion. These mass functions are plotted in Fig. 2. Both 

*It has been suggested that halos identified by spherical overdensity are better suited to observations (although 
which is the most appropriate halo-finder likely depends on the observable) and there is a scatter between FOF and 
spherical overdensity identified halos [102, 96]. However, uncertainties in the mass function from baryonic effects 
appear to be comparable to uncertainties from halo finding [103, 104]. 



13 



the Edgeworth and log-Edgeworth mass functions appear to provide a good fit to simulations in all 
but the most extreme cases; only for /nl = ±500 do we see disagreement. The log-Edgeworth mass 
function appears to do a better job at describing cosmologies with Jnl < ^ and also those with 

In Fig. 3 we show the ratio of the non-Gaussian mass function with ttvl = 2(|/7vl)^ to the one 
with tnl = {^fNhf'-, for a fixed value of Jnl- We find that these two mass functions differ by as 
much as a factor 2 at high redshift and halo mass. Again both the Edgeworth and log-Edgeworth 
mass functions do a reasonable job of predicting the effect of varying tnl^ though the error bars are 
rather large and there seems to be some disagreement towards the high mass end for /tvl = 500. 

Figure 4 shows the non-Gaussian correction to the mass function measured from simulations 
with Jnl = Tnl = and qnl = ±10^ or g^L = ±5 x 10^. As expected, positive qnl increases the 
abundance of massive halos while negative g^L decreases their abundance. For g^i = 10^, both 
the Edgeworth and log-Edgeworth mass functions are in reasonable agreement with simulations, but 
with (77VL = 5 X 10^, the log-Edgeworth mass function is in much better agreement.^ 

As can be seen in Fig. 2 and 4, the improvement from the log-Edgeworth truncation is most 
significant at high masses and redshifts and/or for large non-Gaussianity. For a 3 x 10^^ h~^MQ 
cluster the difference is < 10% at z = even for /jvl = ±500 (and tnl = {(^/^Inl)^)- On the 
other hand, if t^l = '^{'o/^fNif' , the two mass functions differ by 10 — 20% at the same redshift. 
For more modest values of /tvl ( ±100) there is not a significant difference between the two mass 
functions for M < 3 x 10^^ /i~^Mo until z > 1 unless tnl > HInl)'^- 

Finally, in Fig. 5 we compare the effects of /nl, 9nl and tnl on the halo mass function. The 
shape of nNc/nG for the two simulated values of tnl is similar. Indeed, from the analytic mass 
function in Eq. (35) we find that a model with tnl > {^Inl)"^ can be made to look like a model 
with Tnl = (I/a^l)^ at most masses by increasing /nl- On the other hand, g^L and Jnl change 
the shape of the mass function in distinct ways: relative to the /nl models, g^L causes a only slow 
increase/decrease in the abundance of low mass halos but dramatically changes the abundance of 
very massive halos. This suggests that at least in principle, the halo mass function can distinguish 
non-Gaussian initial conditions with skewness and kurtosis (i.e. gNL = with nonzero /nl and tnl) 
from those with kurtosis only (i.e. Jnl = tnl = ^ with nonzero gNL)- This is complementary to 
constraints from large-scale halo clustering, where Jnl and gNL produce roughly degenerate effects 
[35], but Tnl can be constrained independently of /nl by measuring stochasticity in the bias [80, 81]. 

5 Discussion 

In this paper, we have compared semianalytic predictions for the halo mass function to A^-body 
simulations with "generalized local" non-Gaussianity parameterized by Jnl, Qnl and tnl (described 

^There is a caveat: for qnl = 5 x 10®, we find tliat K,fi{M) > K4,{M), wliicfi suggests that tfie 0{gNL) terms in the 
series expansion (34) are not larger than the 0{g%]^) terms which have been neglected. Nevertheless, we find empirically 
that the log-Edgeworth mass function (35) agrees with simulations for these values of qnl- Our perspective is that 
while this point is somewhat unsettling, it is the simulations rather than the Press-Schechter model that validate the 
analytic mass function, so the log-Edgeworth mass function does apply for these qnl values. The current constraints 
on qml are anyway strict enough to avoid this issue. 
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Figure 5: Comparison of the non-Gaussian corrections to the mass function for initial conditions 
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0, 



with fNL = ±500 and tnl = {^Inl) , Jnl = ±500 and tnl 
gNL = ±5 X 10^. The tnl models are rather degenerate; a model with tnl / (f /a^l)^ can be made 
to look like a model with tnl = (f /tvl)^ by adjusting Jnl- On the other hand, g^i and Jnl change 
the shape of nNc/nG in distinct ways. 



in §2). 

Our main result is the log-Edgeworth mass function, Eq. (35), which directly relates the non- 
Gaussian part of the halo mass function to non-Gaussianity in the primordial curvature, by expressing 
{"nNG/f^G) ill terms of the non-Gaussian cumulants K2, K3 and K4 (defined in Eqs. (17) and (20)). 
The log-Edgeworth mass function contains no free parameters and is based on the Press-Schechter 
model for the halo mass function, expanded as a power series in cumulants. A strength of this 
approach, in comparison to a pure fitting function for the non-Gaussian mass function, is that the 
log-Edgeworth mass function can be applied to models of primordial non-Gaussianity not explicitly 
considered in this paper, since the cumulants can be computed from first principles in a given model. 

We have considered the mass function in more detail for two types of non-Gaussian initial con- 
ditions. The first case (§2.1) is a "$^-type" local model in which the coefficients of the three-point 
and four-point functions are independent parameters /atl, tnl- This is implemented by taking the 
initial curvature to be a sum of Gaussian and non-Gaussian fields with constant relative amplitude, 
as discussed recently in [80]. The second case (§2.2) is a "$^-type" model with a ^ATL-term in the 
initial curvature. 

In both of these models, we calculate the cumulants 1^2,1 ^^4 using methods developed in 
Appendix A, and give fitting functions in Eqs. (18)-(21). Thus, for generalized local models with 
parameters Jml-, 9nl, tnl, a completely explicit expression for the non-Gaussian mass function is 
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obtained by plugging these fitting functions for the k's into Eqs. (27)-(32) for the F functions and 
their derivatives, and then into the log-Edgeworth mass function Eq. (35). While we focus on the 
abundance of dark matter halos at z < 2, these expressions may also be useful for other probes 
of non-Gaussianity that are sensitive to the one-point PDF of the density field, e.g. the onset of 
reionization, the abundance of voids or the Lyman-alpha flux PDF [105, 106, 107]. 

In the first model, the cumulants are logarithmically divergent at second order in /tvl (or first 
order in ttvl)- This divergence is a real property of f^L models: if we set up initial conditions via 
$ = + fNLi^G ~ {^g)'^) in a finite box, measurable quantities such as the power spectrum or 
halo mass function will "run" with the box size and slowly diverge in the infinite- volume limit. Our 
interpretation is that /nl models are only defined if a box size is also specified, and our calculation of 
the cumulants includes box size dependence explicitly (which is propagated to box size dependence 
in the log-Edgeworth mass function). 

Keeping track of the box size dependence is potentially a strength of our approach, since the 
statistics of an /tvl model in a Hubble-sized box (the relevant scale for observational constraints) 
may differ from the smaller box sizes typically used to study /nl models using simulations. Using the 
log-Edgeworth mass function, we can predict the ratio n(M)/j^/n(M)^Q of mass functions between 
box sizes Ri = 2 X 10^ h-^ Mpc and Rq = 1600 Mpc. For /nl within the current bounds 
and {tnlI < (few)/'^^, we predict that n{M)ji^ and n(Af)Rg differ by < 1% for halo masses 
< 3 X IO^^/i^^Mq at z = 0. However, even for low values of /tvl we predict an important difference 
if tnl is allowed to be significantly different from {^/nl)"^- For example, if /nl = 20 and tnl ~ 10^ 
we find n{M) / n{M) jif^ ~ 1.1 by M ~ 3x 10^^ and larger at higher masses. We have not attempted 
to compare our predictions for the box size dependence with simulations in this paper, leaving this 
for future work. 

For the first model, we find excellent agreement between the log-Edgeworth mass function and 
A^-body simulations for Jnl G {±250, ±500} and tnl G {{fjNLf, 2(|/ArL)2}. For the second 
model, we find excellent agreement for g^vL G {±1 x 10^, ±5 x 10^} (see Figs. 2-4). These parameter 
ranges are larger than current observational limits and therefore the log-Edgeworth mass function 
appears to provide an excellent fit over the observationally relevant range. 

The log-Edgeworth mass function constructed in this paper is conceptually similar to the Edge- 
worth mass function from [46]; the main difference is that we expand the quantity ln(F(M)) defined 
in Eq. (22) as a power series in cumulants, rather than the mass function n(M). The log-Edgeworth 
and Edgeworth mass functions agree in the limit where many terms in the series expansion are 
retained (or the limit of small /nl) but differ in practice if the expansions are truncated at a fixed 
finite order. The main difference is that the Edgeworth mass function has non-physical asymptotic 
behavior at large M, but in the observationally relevant range the log-Edgeworth mass function con- 
tinues to behave correctly in this limit. Additionally, even away from the large-M limit, we find that 
the log-Edgeworth mass function is a slightly better fit to the simulations than the Edgeworth mass 
function. Recent hints of an overabundance of very massive clusters have renewed interest in un- 
derstanding the effects of primordial non-Gaussianity on the abundance of rare objects. Accurately 
predicting the high-mass, high redshift end of the halo mass function is critical to interpreting the 
significance of systems. The log-Edgeworth mass function appears to describe the halo abundance 
from simulations in this regime. 
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It would be interesting to compare the framework developed in this paper with simulations in 
more general models, particularly "single-field" non-Gaussianity with parameters /^^'^ and f^Jj^ 
[22]. The smoothed cumulants of the density field don't retain information about the "shape" 
of the A^-point functions and one may therefore expect this framework to work for more general 
non-Gaussian initial conditions. (The fitting functions for the kat's would of course need to be 
recomputed.) Since the log-Edgeworth mass function correctly predicts the mass function over 
a wide range of the {fNL,gNL,TNL} parameter space considered in this paper, one has increased 
confidence that it may be more applicable to more general non-Gaussian initial conditions. However, 
we defer a detailed comparison to future work. 
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A Calculating cumulants 

In §3.1, we quoted fitting functions for the following cumulants: the ©(/atl) contribution to AC3, 
the 0{gNL) contribution to K4, and the 0{tnl) contributions to K2 and K4. In this appendix, we 
describe our procedure for obtaining these results, via a combination of analytic and Monte Carlo 
methods. 

A.l Convergent cumulants: K^{M)fj^^ and K4(M)g^^ 

We first calculate K^{M)fj^^ and K4^{M)g^^. These are the cumulants which do not diverge in 
the infinite- volume limit. In these cases, we can compute the cumulant efficiently by numerical 
integration, as we now explain. 

The smoothed linear overdensity field 6m and the curvature $ are related by 



where we have introduced the notation aM{k) = W{kR{M))a{k). (We have ignored redshift 
dependence, since it will eventually drop out when we compute reduced cumulants of the form 



{^M{^f)fML = ^fNL / j7^j7^j7^aM{ki)aM{k2)aMik3)P,s>{ki)P^{k2){27Tf6D{^^ 



The 6D integral can be rewritten in a more tractable form as follows. We write the delta function 

^"Note that we use a notation in which contributions of different orders are distinguished by subscripts, e.g. the first 
hne of Eq. (f9) is denoted K4(Af)g„j^ and the second hue is denoted /s;4(Af)Tjvi- 




(37) 




(38) 
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as an integral (27r)'^(^£)(^ kj) = J d?Te ^E^^*'^ and introduce the notation: 



aM{r) = I aM{k)e 



(27r)3 

A:2aAf(fc)jo(fcr) (39) 







27r2 







^k^aM{k)P^{k)jo{kr) (40) 



where jo{kr) = sm{kr)/{kr) and we have used spherical symmetry to write each Fourier transform 
as a ID integral. Then 



oo 



{5M{^Y)h^ = ^fNL / drA7Tr'dM{r)(3M{ry (41) 

Jo 

and (t(M)2 = J dr ATTr'^aM{f)/3M{f), so that the reduced cumulant is given by: 

dr ATTr'^aM{r)l3M{rf 
47rr2aM (?^)/3m (?•))' 

and in this form, the cumulant can be computed efficiently, as a sequence of ID integrals. We note 
that this procedure can be extended to find an expression for the /jvl ^ contribution to in terms 
of a sequence of ID integrals. 

A similar calculation can be done for the cumulant K4(Af)gjvL- 

2457VL /" d^ki d^\L2 d^kg d^lii 
^(Mj4 J (27r)3 (27r)3 {2^f {2tt) 

xP^{h)P^{k2)P^ik3){27Tf6D{lil + k2 + kg + k4) 



'^'i{M)fNL = ^fNL , _ , (42) 



'^4(M)3^^ = — — ^ / .^^.3 .^^.3 .3 ,^.^ aM{ki)aM{k2)aM{kz)aM{ki) 
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/o°°rfr47rr2dAf(r)/3M(r) 



47rr2dA/(r)/3A/(?')) 
A. 2 Infrared-divergent cumulants: H2{M)tnl K4{M) 



The cumulants K2(M)t-^^ and K4(M)t-jv^ are formally divergent in the limit where the simulation 
volume goes to infinity, or equivalently as integrals over wavenumbers kj are computed with lower 
cutoffs /cmin " ^ 0. This divergence can be understood physically as follows (for more discussion see 
[108]). As fcmin — )• 0, new long-wavelength modes of "appear", which can couple to modes at a 
fixed physical scale k (via the term in Eq. (2)) and generate divergent contributions to *^(k). 
For example, the power spectrum P^{k) contains a log-divergent term proportional to /^^. This 
divergence can be regulated by fixing a box with finite volume and requiring that = Oi 

where the mean (•) is defined by integrating over the box. For a fixed L, all cumulants will be finite, 
but some contributions diverge as ln(L) as L — ?• oo. 

Our perspective is that /nl models are only defined if a box size L is also specified. When 
comparing with A^-body simulations, we will choose L to be the size of the simulation box. This 
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choice is consistent with the way the A^-body initial conditions are generated: we do not simulate 
power in the DC mode of the simulation box, and this corresponds to regulating the divergence 
using a volume L^. When comparing to observations, we suggest choosing L to be large enough to 
enclose the Hubble volume. 

The correction to the variance K2(M)t-^^ is given by 

-2(MV^, = / + ^^2\fP^{h)P^{h) (44) 

and this integral diverges in the limits /ci — t- and A;2 — 0. In a finite box of volume L^, it can 
be regulated by replacing / — )■ -p- X^k^c where the sum ranges over discrete Fourier modes 
(i.e. modes of the form {kx, ky, k^) = {2Trnx/L, 2TTny/L, l-Kriz/L) where the rii are integers). 
Next we consider the ttvl contribution to the kurtosis, 



TNL 1 f d^^i d^k2 d^kg d^ki 

(6/5)2 ^(^/)4 J (2^(2^(2^(2^ 

xP^{\ki + k2\)P^{ki)P^{ks){27rf 60(^1 + k2 + k3 + k4) , (45) 



K4(M)rjvi = 48 / v^ QAf (fcl)aAf(fc2)QA/(fc3)«Af (^4 



this integral contains a divergence as the internal momentum q = (ki + k2) approaches zero. To see 
this explicitly, we change variables: 



(6/5)2 ^(j^)4 J (27r 



d^k 

(2^ 



QM(A:)aA/(|q - k|)P$(/c) 



(46) 



in the q — )■ limit, the expression in brackets approaches it(M)2 = J -^^aM{k)'^P^ik), a nonzero 

(and finite!) value. The divergence can be regulated by replacing f ^^^^3 — L^J2ci:^o- 

It is useful to calculate the infrared divergent part of these cumulants, i.e. the leading behavior 
in the infinite-volume limit. In both cases, this will follow from analyzing the infared divergent part 
of the integral / -^^P^{k), after regulating by replacing the integral by a discrete sum over Fourier 
modes in a finite volume L^. For a scale-invariant power spectrum of the form P^{k) = 27r2A|/A;^, 
it is easy to see that 

/d^k 
j^P'S'ik) = A| ln(L) + (finite) . (47) 

Comparing with Eqs. (44) and (46), we get the IR-divergent terms: 



K2(M).^, = 4^^A2ln(L) + (finite) (48) 



«4(M).^, = 48^^A2ln(L) + (finite). (49) 

Note that the IR-divergent part is independent of M in both cases. 

The preceding expressions assume scale invariance for simplicity but if < 1 , we find empirically 
that an excellent fitting function for the IR-divergent behavior is given by: 

which agrees with Eq. (47) in the limit — )• 1. 
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A. 3 Monte Carlo simulations 



We now describe a method for efficiently estimating cumulants for a fixed box size, via Monte Carlo 
simulations of the density field. 

In each Monte Carlo realization, we simulate a Gaussian initial curvature define fields 

^M, SIj, 61} by 

<5Af(k) = aM{k) I d^xe-^'^-^'^Gi^) 

6Uk) = aM{k) J d3xe"^i^ ''cDc(x)2 

= aM{k)ld'^e~^''-^{^G{^f-3{^l)^G{^)). (51) 

The field 6m represents the linear density field smoothed on mass scale M, and the fields 6^^, 
represent non-Gaussian contributions of Jnl-^YP^ or ^ArL-type respectively. 
We then estimate cumulants by 

.3(M),,, = 3f^, (5,,(x)2)3/2 
(M\ A {hl{^f6rj{^)) 



i2\ 



AJ2(M) 



TNL {6*M{yi 

(6/5)2 (^^^(x)2 



-4(M).,, = ^6/5)2 (^,Kx)^)^ ^''^ 

where (•) denotes an average over MC realizations and also over position x within each realization. 

This Monte Carlo scheme isolates contributions of a given order (e.g. the order-jAr^, contribution 
to K^{M) is estimated without any contribution from the order- term) and results in good 
computational efficiency. For example, one can get reasonable-looking results after a single Monte 
Carlo simulation. 



A. 4 Fitting functions 

The convergent cumulants ks{M) and Ki{M)gj^^ can be calculated either by numerical integration 
(§A.l) or by Monte Carlo simulations of the density field (§A.3). We checked that the two agree, 
and find that the fitting functions given previously in Eq. (18) and the first line of Eq. (19) are 
excellent approximations (see Fig. 1). 

The divergent cumulants K2{M)t-^^ and ^^{M)^-^^^ are more subtle since dependence on the box 
size L must be included in the fitting functions. In these cases, the IR-divergent pieces can be 
calculated in closed form (§A.2), but calculating the precise values of the cumulants at a finite value 
of L is awkward since the integrals must be regulated by summing over a discrete set of Fourier 
modes. We obtain fitting functions using a hybrid approach: we use the Monte Carlo scheme from 
§A.3 with fixed box size Lq = 1600 h^^ Mpc, and find an empirical fitting function for the M 
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dependence. We then add a term proportional to ln(L/Lo) with prefactor chosen to match the 
IR-divergent term calculated in Eq. (49). The fitting functions for K2(^)tjvl ^^'^ ^4^i^)TML given 
previously in Eq. (21) and the second line of Eq. (19) were obtained by this procedure. As a check, 
we find that these fitting functions are excellent approximations for box sizes Lo/4 < L < Lq. Since 
they have been constructed in a way which guarantees correct asymptotic behavior as L — )■ oo, 
we anticipate that they can be safely extrapolated to L ^ Lq, where doing simulations would be 
computationally prohibitive. This property is important since L should be chosen to be Hubble-sized 
when comparing with observations, as discussed in §3.1. 
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